---
title: "Replication Code"
subtitle: "Who supports subsidizing renewable energy via an electricity surcharge? Evidence from an online survey in South Korea"
date: "`r Sys.Date()`"
output: pdf_document
editor_options: 
  chunk_output_type: console
---

# 1. Preparation

```{r}
pacman::p_load(MASS, readxl, tidyverse, car, rethinking, coda, simcf, tile)
source("theme_caviz.R")

data <- read.csv("data_1123.csv")[,-1]

data$pop <- as.numeric(data$pop)
data$popperkm2 <- data$pop / data$area
data$capaperha <- 100 * data$capa_sum / data$area

data$q1 <- recode(data$q1,"1=5;2=4;4=2;5=1")
data$q2 <- recode(data$q2,"1=5;2=4;4=2;5=1")
data$q3 <- recode(data$q3,"1=5;2=4;4=2;5=1")
data$q4 <- recode(data$q4,"1=5;2=4;4=2;5=1")
data$q5 <- recode(data$q5,"1=5;2=4;4=2;5=1")
data$q6 <- recode(data$q6,"1=5;2=4;4=2;5=1")

varnames <- c("id","support","impact","pollute","reliable",
              "howfreq","landscape","income","educ",
              "ideology")

colnames(data)[1:10] <- varnames

data$logpopden <- round(log(data$popperkm2+1),3)
data$logcapden <- round(log(data$capaperha+1),3)
```

# 2. Descriptive Statistics

```{r}
# Table 1
mean(data$age);sd(data$age)
mean(data$income);sd(data$income)
mean(data$educ);sd(data$educ)
mean(data$ideology);sd(data$ideology)
data %>% group_by(gender) %>% count() %>% mutate(percent=100*n/1043)
data %>% group_by(income) %>% count() %>% mutate(percent=100*n/1043)
data %>% group_by(educ) %>% count() %>% mutate(percent=100*n/1043)
data %>% group_by(ideology) %>% count() %>% mutate(percent=100*n/1043)
data %>% group_by(support) %>% count() %>% mutate(percent=100*n/1043)
```

```{r}
# Figure prep
data %>% group_by(support) %>% count() %>% mutate(percent=100*n/1043) %>% 
  ggplot(aes(x=support, y=percent)) +
  geom_bar(stat="identity") +
  xlab("How much do you agree with increasing the renewable energy surcharge?") +
  ylab("Responses (%)") +
  scale_x_discrete(limit=c(1,2,3,4,5),
                   labels=c("Strongly disagree","Disagree",
                            "Neither disagree \n nor agree",
                            "Agree","Strongly agree")) +
  theme_caviz_hgrid +
  theme(aspect.ratio = 1/1.1618, 
        axis.text.x = element_text(color = "black", size = 11),
        axis.ticks.x = element_line(size=1),
        axis.title.x.bottom = element_text(size=15),
        axis.title.y = element_text(size=13))

# Figure 2

data %>% group_by(support, landscape) %>% count() %>% mutate(percent=100*n/1043) %>% 
  ggplot(aes(x=support, y=percent, fill=factor(landscape))) +
  geom_bar(stat="identity") +
  xlab("How much do you support/oppose with increasing the renewable energy surcharge?") +
  ylab("Responses (%)") +
  scale_x_discrete(limit=c(1,2,3,4,5),
                   labels=c("Strongly oppose","Oppose",
                            "Neither oppose \n nor support",
                            "Support","Strongly support")) +
  labs(fill="Do you agree that solar farms disrupt surrounding landscapes?") +
  scale_fill_manual(
    labels=c("Strongly disagree","Disagree","ND/NA","Agree","Strongly agree"),
    values=c("1"="#5A95EA", "2"="#9FBDE8", "3"="#FAF4C9","4"="#F0BDA7","5"="#EB6B34"))+
  theme_caviz_hgrid +
  theme(legend.position = "top",
        legend.title = element_text(size=11),
        legend.text = element_text(size=11),
        aspect.ratio = 1/1.1618, 
        axis.text.x = element_text(color = "black", size = 11),
        axis.ticks.x = element_line(size=1),
        axis.title.x.bottom = element_text(size=15),
        axis.title.y = element_text(size=13))

# Figure 3

data %>% group_by(support, pollute) %>% count() %>% mutate(percent=100*n/1043) %>% 
  ggplot(aes(x=support, y=percent, fill=factor(pollute))) +
  geom_bar(stat="identity") +
  xlab("How much do you agree with increasing the renewable energy surcharge?") +
  xlab("How much do you support/oppose with increasing the renewable energy surcharge?") +
  ylab("Responses (%)") +
  scale_x_discrete(limit=c(1,2,3,4,5),
                   labels=c("Strongly oppose","Oppose",
                            "Neither oppose \n nor support",
                            "Support","Strongly support")) +
  labs(fill="Do you agree that solar energy produces less air pollution than fossil fuel?") +
  scale_fill_manual(
    labels=c("Very disagree","Disagree","ND/NA","Agree","Very agree"),
    values=c("1"="#B249EE", "2"="#CB93EA", "3"="#FAF4C9","4"="#8DBFA1","5"="#199149"))+
  theme_caviz_hgrid +
  theme(legend.position = "top",
        legend.title = element_text(size=11),
        legend.text = element_text(size=11),
        aspect.ratio = 1/1.1618, 
        axis.text.x = element_text(color = "black", size = 11),
        axis.ticks.x = element_line(size=1),
        axis.title.x.bottom = element_text(size=15),
        axis.title.y = element_text(size=13))

# Figure 4

data %>% group_by(support, impact) %>% count() %>% mutate(percent=100*n/1043) %>% 
  ggplot(aes(x=support, y=percent, fill=factor(impact))) +
  geom_bar(stat="identity") +
  xlab("How much do you agree with increasing the renewable energy surcharge?") +
  xlab("How much do you support/oppose with increasing the renewable energy surcharge?") +
  ylab("Responses (%)") +
  scale_x_discrete(limit=c(1,2,3,4,5),
                   labels=c("Strongly oppose","Oppose",
                            "Neither oppose \n nor support",
                            "Support","Strongly support")) +
  labs(fill="Do you agree that climate change is negatively affecting your daily lives?") +
  scale_fill_manual(
    labels=c("Very disagree","Disagree","ND/NA","Agree","Very agree"),
    values=c("1"="#951222", "2"="#CF6671", "3"="#FAF4C9","4"="#6898DA","5"="#0C54B8"))+
  theme_caviz_hgrid +
  theme(legend.position = "top",
        legend.title = element_text(size=11),
        legend.text = element_text(size=11),
        aspect.ratio = 1/1.1618, 
        axis.text.x = element_text(color = "black", size = 11),
        axis.ticks.x = element_line(size=1),
        axis.title.x.bottom = element_text(size=15),
        axis.title.y = element_text(size=13))


```


# 3. Bayesian Ordinal Model (baseline)

```{r}
set.seed(2022)
bayesdata <- list(
  S=data$support, # dependent variable
  L=data$landscape, # landscape variable
  P=data$pollute, # pollution variable
  I=data$impact, # impacts from climate variable
  income=data$income, educ=data$educ, ideology=data$ideology, 
  age=data$age, gender=data$gender,
  capaden=data$logcapden, popden=data$logpopden, forest=data$forestper,
  alpha=rep(2,4)
)
```


$$
\begin{aligned}
\textrm{Support}_i & \sim \textrm{Ordered-logit}(\phi_{i}, \kappa) \\
\phi_{i} & = \beta_{L} \sum_{j=0}^{L_i-1} \delta_{Lj} +  \beta_{P} \sum_{j=0}^{P_i-1} \delta_{Pj} + \beta_{I} \sum_{j=0}^{I_i-1} \delta_{Ij}  \\
& + \beta_1 \textrm{income}_i + \beta_2 \textrm{educ}_i  + \beta_3 \textrm{ideology}_i  + \beta_4 \textrm{age}_i  + \beta_5 \textrm{gender}\\
\kappa_k & \sim \textrm{Normal}(0, 1.5) \\
\beta_L,\beta_P,\beta_I,\beta_1,\dots,\beta_6 & \sim \textrm{Normal}(0,1) \\
\delta_{Lj},\delta_{Pj},\delta_{Ij} & \sim \textrm{Dirichlet}(2,2,2,2)
\end{aligned}
$$

```{r}
model1 <- ulam(
  alist(
    S ~ ordered_logistic(phi, kappa),
    phi <- bL*sum(delta_L[1:L]) +  bP*sum(delta_P[1:P]) + bI*sum(delta_I[1:I]) +
      b1*income + b2*educ + b3*ideology + b4*age + b5*gender,
    kappa ~ normal(0, 1.5),
    c(bL,bP,bI,b1,b2,b3,b4,b5) ~ normal(0, 1),
    vector[5]: delta_L <<- append_row(0,deltaL),
    vector[5]: delta_P <<- append_row(0,deltaP),
    vector[5]: delta_I <<- append_row(0,deltaI),
    simplex[4]: deltaL ~ dirichlet(alpha),
    simplex[4]: deltaP ~ dirichlet(alpha),
    simplex[4]: deltaI ~ dirichlet(alpha)
  ), data=bayesdata, chains=4, cores=4, iter=5000
)
```

```{r}
precis(model1, depth=2, digits=3, omit="kappa", prob=0.95)

# Figure 5
model1_result <- precis(model1, depth=2, digits=3, omit="kappa", prob=0.95)
coefs <- model1_result[c(8,7,6,5,4,3,2,1),]
coefs <- round(coefs, digits=2)
rownames(coefs) <- c(
  "Landscape disruption", "Air quality", "Impacts of climate change",
  "Household income", "Education", "Ideology", "Age", "Gender")

labels <- c("Landscape disruption \n [-1.20, .00]",
            "Air quality \n [1.27, 2.51]",
            "Climate change impact \n [2.68, 4.31]",
            "Household income \n [-.02, .14]", 
            "Education \n [-.08, .14]", 
            "Ideology \n [-.48, -.17]", 
            "Age \n [-.02, .00]", 
            "Gender \n [-.24, .21]")

coef_plot <- ropeladder(x=coefs$mean,
                        lower=coefs$`2.5%`,
                        upper=coefs$`97.5%`,
                        labels=labels, plot=1, mark="dashed", size=0.5)
tile(coef_plot,
     limits=c(-1.5, 4.5),
     gridlines=list(type="xt", add=T, lty="solid"),
     topaxis=list(add=T, at=seq(-1.0, 4.0, by=1.0),
                  labels=c("-1.0","0","1.0","2.0","3.0","4.0")),
     xaxis=list(add=T, at=seq(-1.0, 4.0, by=1.0),
                  labels=c("-1.0","0","1.0","2.0","3.0","4.0")),
     xaxistitle=list(labels="Parameter value"),
     width=list(null=4),
     height=list(xaxistitle=3, plottitle=1),
     output=list(outfile="figure5", width=6, height=6/1.618))

```

```{r}
pairs(model1, pars="deltaP")
```


# 4. Robustness Check (including geographical factors)

$$
\begin{aligned}
\textrm{Support}_i & \sim \textrm{Ordered-logit}(\phi_{i}, \kappa) \\
\phi_{i} & = \beta_{L} \sum_{j=0}^{L_i-1} \delta_{Lj} +  \beta_{P} \sum_{j=0}^{P_i-1} \delta_{Pj} + \beta_{I} \sum_{j=0}^{I_i-1} \delta_{Ij}  \\
& + \beta_1 \textrm{income}_i + \beta_2 \textrm{educ}_i  + \beta_3 \textrm{ideology}_i  + \beta_4 \textrm{age}_i  + \beta_5 \textrm{gender}_i \\
& + \beta_7 \textrm{capaden}_i + \beta_8 \textrm{popden}_i + \beta_9 \textrm{forest}_i \\
\kappa_k & \sim \textrm{Normal}(0, 1.5) \\
\beta_L,\beta_P,\beta_I,\beta_1,\dots,\beta_9 & \sim \textrm{Normal}(0,1) \\
\delta_{Lj},\delta_{Pj},\delta_{Ij} & \sim \textrm{Dirichlet}(2,2,2,2)
\end{aligned}
$$

```{r}
model2 <- ulam(
  alist(
    S ~ ordered_logistic(phi, kappa),
    phi <- bL*sum(delta_L[1:L]) +  bP*sum(delta_P[1:P]) + bI*sum(delta_I[1:I]) +
    + b1*income + b2*educ + b3*ideology + b4*age + b5*gender 
    + b6*capaden + b7*popden + b8*forest,
    kappa ~ normal(0, 1.5),
    c(bL,bI,bP,b1,b2,b3,b4,b5,b6,b7,b8) ~ normal(0, 1),
    vector[5]: delta_L <<- append_row(0,deltaL),
    vector[5]: delta_P <<- append_row(0,deltaP),
    vector[5]: delta_I <<- append_row(0,deltaI),
    simplex[4]: deltaL ~ dirichlet(alpha),
    simplex[4]: deltaP ~ dirichlet(alpha),
    simplex[4]: deltaI ~ dirichlet(alpha)
  ), data=bayesdata, chains=4, cores=4, iter=5000
)
```

```{r}
precis(model2, depth=2, digits=3, omit="kappa", prob=0.95)

model2_result <- precis(model2, depth=2, digits=3, omit="kappa", prob=0.95)
coefs2 <- model2_result[c(11,9,10,8,7,6,5,4,3,2,1),]
coefs2 <- round(coefs2, digits=2)
rownames(coefs2) <- c(
  "Landscape disruption", "Air quality", "Impacts of climate change",
  "Household income", "Education", "Ideology", "Age", "Gender",
  "Solar farm density", "Population density", "Forest coverage")
```


